13. Zero-Inflated Beta Models for Microbiome Data

13.2.1 Introduction to ZIBseq

ZIBSeq (Peng et al. 2016) was developed to identify differentially abundant features between multiple clinical conditions while addressing the unique characteristics of metagenomics data with small sample size, high dimensionality, sparsity, often with a large number of zeros and skewed distribution under the compositions (proportions) setting. ZIBSeq directly handles the proportion data, assuming that proportional dependent variable can be characterized by the beta distribution.

13.2.1.1 ZIBseq Model

Assume x follows a beta distribution: x~Beta(μ, ϕ), then the density of beta distribution (Ferrari and Cribari-Neto 2004) can be described as a function of u and ϕ:

(13.1)

The mean and variance can be parameterized as:

(13.2)

where Γ(⋅) is the gamma function, μ(0 < μ < 1) is the mean, and ϕ(ϕ > 0) is a precision parameter.

Beta distribution has a wide range of shapes depending on the values of mean and precision parameters. Thus, beta regression models are very suitable to model the continuous response variables of rates and proportions when their values are restricted to the interval (0,1). To model the variable containing zero or one, Ospina and Ferrari in 2012 (Ospina and Ferrari 2012) proposed a more general class of zero-or-one inflated beta regression models for continuous proportions. Given that microbiome data often contains zero instead of one, ZIBSeq was developed based on the zero-inflated beta regression, assuming the response variable has a mixed continuous-discrete distribution with probability mass at zero.

The zero-inflated beta distribution adds a new parameter α to account for the probability of observations at zero. The subsequent mixture density is given as below:

(13.3)

where f(x; μ, ϕ) is the beta density in (13.1). The zero-inflated beta and hence ZIBSeq directly model the proportion of the response variables. Thus, the feature abundances are first normalized into the proportion of features. Let denote the normalized feature abundance, i.e., the proportion of feature j reads in sample i, ZIBSeq simply normalizes the feature abundance by dividing read counts of each given feature by the total mapped read counts in the sample: . Thus, the obtained are proportions of feature j, j = 1, …, m on n samples. By assuming that each has a probability density function (13.3) with parameters , , and . To fit the parameters in mixture distribution (13.3) for feature j, the ZIBSeq model is defined as below:

(13.4)

In the above equations, , , and are unknown regression parameters to be estimated, while y_i is an outcome measurement indicating the class of sample i, while T_i is the ith sample depth. The dispersion parameter ϕ in (13.4) is approximated by a binomial distribution and a beta distribution under parameterization in (13.1). First, assume the number of reads on a particular feature c follows a binomial distribution with the sample depth T and the probability that the sample reads in this feature μ. Second, let , then, the variance of proportion x will be . Third, in the case, the proportion x follows a beta distribution under parameterization in (13.1), the variance will be . An approximation of the dispersion parameter ϕ: ϕ ≈ T − 1 is derived from these two forms of variance.

13.2.2 Implement ZIBseq

The microbiome data structure that ZIBseq works on is a n × m matrix (see Table 11.​1), with each of n rows presenting samples collected from two or more classes, and each of m columns presenting expression levels of genes or abundance of taxa measured for each sample. Each entry c_ij denotes the number of reads from sample i that mapped to feature (i.e., taxon) j. denotes the total number of reads of sample i. Each count c_ij not only depends on the reads of taxon j but also varies from the sequencing depth: the total number of reads of sample i. In Table 13.1, y_i, i = 1, …, n denotes the outcome variable associated with sample i, and it could be class membership of the sample.

Table 13.1

Microbiome data structure for ZIBseq model

	Feature1	Feature2	…	Featurem	Total	Outcome
Sample1	c11	c12	…	c1m	T1	y1
Sample2	c21	c22	…	c2m	T2	y2
…			…	…	…	…
Samplen	cn1	cn2	…	cnm	Tn	yn

After running this function, the significant feature (“sigFeature”), the features that have been concerned (“useFeature”), and both Q-value (“qvalue”) and P-value (“pvalue”) will be provided in the output.

The “otu_table_fecal_genus_vdr.csv” is a csv file containing 16S rRNA sequencing data matrix for the study samples. The VDR mice abundance data provides the absolute sequence read counts for a total of 240 taxa at the genus rank level (rows) in 8 samples (columns), preprocessed as described in the main article. Among these eight samples, five samples with VDR were knockout and the remaining three samples were wild types. We are interested in whether dysbiosis and functions of intestinal microbiome are associated with the VDR knockout. Thus, the main metadata is the VDR condition status. We can provide additional metadata, which could be a CSV file containing the sample metadata in the data matrix to run the ZIBSeq. Here, we use commands to extract the grouping information from the sequencing data matrix.

The data frame has 249 columns, the first column is grouping variable, the remaining 248 columns are taxa.

Of course, we can create a meta dataset a priori and then load it into R and merge it to taxa abundance dataset.

Then, we can perform ZIBSeq analysis. The following steps including screening features, normalizing and transforming data as well as correcting multiple hypothesis testing are processed automatically by the ZIBSeq() function or through specifying such as to specify transform = TRUE or T will do a square root or cube root transformation. We explicitly list them below for reminding us what are really the function ZIBSeq() working on.

Because taxa with a small number of reads are not reliable, so the function first remove any features with total counts less than 2 time large of the sample size to reduce the effects of noises and measurement errors.

ZIBSeq uses a simple normalization procedure to convert the raw abundance measure to a proportion by dividing each feature read count by the total feature read counts in the sample, which results in relative abundance measure ranging [0,1). After normalization, a square root or cube root transformations is performed to ensure the proportion data are better fitting a beta distribution if the distributions of the proportion are extremely left skewed.

13.2.3 Remarks on ZIBSeq

13.3.1 Introduction to ZIBB

ZIBB (Hu et al. 2018) was proposed based on zero-inflated beta-binomial model to detect the association of microbiome with a continuous or categorical phenotype of interest. ZIBB model assumes that the count data distribution consists of a mixture of a point mass probability at zero (zero model) and a beta-binomial distribution (count model). Thus, ZIBB model has two components: (1) a zero model accounting for excess zeros and (2) a count model (the remaining component) allowing for over-dispersion effects through beta-binomial regression, which may have an appreciable additional mass at zero. The proposed ZIBB method is an extension of the beta-binomial model of BBSeq (Zhou et al. 2011) that was developed for the analysis of RNA sequence count data. Employing zero-inflated modeling and considering a constrained approach to estimate the over-dispersion parameters have been considered as two major improvements in the proposed ZIBB framework (Hu et al. 2018).

13.3.1.1 Zero and Count Models in ZIBB

ZIBB was proposed for analysis of 16S rRNA sequence count data, a m × n sequence count matrix with columns representing samples and rows representing OTUs (ASVs/Taxa). Let Y = (y_ij) ∈ Ζ^m × n be the count matrix, and each element y_ij represents the count of OTU i in sample j, where (i = 1, …, m, j = 1, …n). Let be the library size for sample j. Then, the distribution of microbiome counts data is modeled by a two-stage canonical beta-binomial model, which assuming y_ij follows a binomial distribution Bin(s_j, μ_ij), and μ_ij~Beta(α_1ij, α_2ij). The probability mass function of the beta-binomial distribution with parameter (s_j, α_1ij, α_2ij) is written as f(⋅| α_1ij, α_2ij) or

(13.5)

In the above formulation, the probability is impliedly dependent on s_j. Since μ_ij~Beta(α_1ij, α_2ij), then the mean and variance of Beta(α_1ij, α_2ij) distribution can be written as:

• E(μ_ij) = α_1ij/(α_1ij, α_2ij) and ϕ_iE(μ_ij)(1 − E(μ_ij)). And it is trivial to solve that α_1ij = E(μ_ij)(1 − ϕ_i)/ϕ_i and α_2ij = (1 − E(μ_ij))(1 − ϕ_i)/ϕ_i.

Through reparameterizing the parameters (α_1ij, α_2ij) of the beta distribution to (E(μ_ij), ϕ_i), ϕ_i ≥ 0 in the reparameterized form can be clearly interpreted as the over-dispersion parameter. ZIBB model is a mixture of a point mass at zero (zero model) and a beta-binomial distribution (count model). The density of count y_ij is written as:

(13.6)

where π_ij is the point mass at zero and f(⋅) is the probability mass function of the canonical beta-binomial distribution as in (13.5). In the above ZIBB model, zero inflation has been taken accounted for by π_ij. Two link functions are used to include the effects of phenotype and covariates in the modeling: one for the zero model and another for the count model. The link function of the zero model is given below:

(13.7)

where z = (z₁, …, z_n)^T ∈ R^n × q is referred to as the design matrix for the zero model. z = (z_{0, j}, …, z_{q − 1, j})^T ∈ R^q is the vector of zero-inflation related covariates for sample j (the z_{0, j} = 1 denotes the intercept) and η_i = (η_{0, i}, …, η_{q − 1, i})^T ∈ R^q is the vector of corresponding coefficients for OTU i. When q = 2, the zero inflation-related covariates are . The link function of the count model is given below:

(13.8)

where the X = (X₁, …X_n)^T ∈ R^n × p is referred to as the design matrix for the count model. X_j = (X_{0, j}, …, X_{p − 1, j})^T ∈ R^p is the vector of phenotypes of interest (X_{0, j} = 1 is the intercept of design matrix) for sample j and β_i = (β_{0, i}, …, β_{p − 1, i})^T ∈ R^p is the vector of corresponding coefficients for OTU i.

To model the relationship between the mean and the variance, as well as between the mean and the over-dispersion. ZIBB uses the polynomial fit:

(13.9)

as the constraint between over-dispersion parameter ϕ_i and coefficients β_i in the count model.

The strategy of considering a constrained approach to model the over-dispersion as a polynomial function of the systematic (mean) component of the generalized linear model was adopted from BBSeq (Zhou et al. 2011). Based on the authors of ZIBB, K = 3 is typically sufficient to effectively model the relationship.

The distinctive feature of ZIBB lies on its modeling the mean-over-dispersion relationship. Thus, the authors of ZIBB called their model and the associated estimation approach as the constrained model and name the model that estimates ϕ_i separately for each OTU i as the free model.

13.3.2 Implement ZIBB

ZIBB for the discovery analysis is implemented via the R package ZIBBSeqDiscovery. ZIBBSeqDiscovery uses zero-inflated beta-binomial model to analyze microbiome count data to detect the associate between phenotype of interest and the composition of the counts. The interested covariates can be adjusted via the link functions of zero and count models. The moment corrected correlation (MCC) approach is applied to adjust the P-values.

ZIBBSeqDiscovery takes two approaches to fit the ZIBB model: (1) The free approach treats the over-dispersion parameters for OTUs as independent, and (2) the constrained approach uses a mean-over-dispersion relationship to the count data.

In the ZIBBSeqDiscovery package, the covariate of interest is referred to as phenotype. Current version (latest version 1.0, March 2018) of this package only considers one phenotype. In the design matrix X, the first column is the intercept, and the phenotype of interest is assumed to corresponding to the second column.

Since the authors of ZIBB (Hu et al. 2018) showed that constrained modeling of ZIBB is better than its free modeling of ZIBB in terms of the type I errors and power, thus we here only report the constrained version of ZIBB using this dataset.

The ZIBB methods may fail to converge for a small proportion of OTUs. This typically occurs when OTUs have few (e.g., four or fewer) non-zero counts and for such OTUs. In such case, the power to detect association with experimental variables is small; thus before modeling using ZIBB methods, those OTUs need to be filtered out. For the remaining small number of OTUs that fail to converge (usually 2% or fewer), the ZIBB methods substitute P-values from the MCC package.

In below analysis, we removed all OTUs with zero counts across all samples. Based on the practice, for the discrete case, OTUs with zero counts across any one of groups were also removed.

We use below commands to construct the design matrix and zero inflation related matrix.

NA mostly appears when the OTU counts are zero in most of the samples. In this case, there are 185 samples, thus we want to check the cases such that OTU has 180, 181, 182, 183, and 184 zero counts across the 185 samples, respectively. The following commands can be used to check the effects of MCC adjustment. The NA counts in the reported P-values before and after MCC adjustment will be printed out.

13.3.3 Remarks on ZIBB

Overall, this study (Hu et al. 2018) showed that both zero inflation and the proposed constraint approaches are vital for accurate and powerful analysis of microbiome count data and ZIBB modeling with the constrained approach is preferred among the competing methods, including BBSeq, edgeR, and ZINB. Especially, it was demonstrated (Hu et al. 2018) that in terms of Type I errors and power for both discrete and continuous cases: (1) ZIBB has higher or comparable performance compared to competing methods. (2) Both the approaches of free and the constrained ZIBB modeling can provide accurate standard error estimates. And (3) the constrained modeling of ZIBB outperformed free modeling of ZIBB. However, the ZIBB methods may fail to converge for a small proportion of OTUs. This converge problem typically occurs when OTUs have four or fewer non-zero counts (Hu et al. 2018).